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Abstract. We study numerically the semiclassical limit for the nonlinear 
Schrodinger equation thanks to a modification of the Madelung transform due 
to E. Grenier. This approach is naturally asymptotic preserving, and allows 
for the presence of vacuum. Even if the mesh size and the time step do not 
depend on the Planck constant, we recover the position and current densities 
in the semiclassical limit, with a numerical rate of convergence in accordance 
with the theoretical results, before shocks appear in the limiting Euler equa- 
tion. By using simple projections, the mass and the momentum of the solution 
are well preserved by the numerical scheme, while the variation of the energy is 
not negligible numerically. Experiments suggest that beyond the critical time 
for the Euler equation, Grenier's approach yields smooth but highly oscillatory 
terms. 



1. Introduction 
We consider the cubic nonlinear equation 

e 2 

(1.1) ied t u e + —Au £ = \u £ \ 2 u £ , (i,x)eR+xR d . 

The goal is to compute the solution u e in such a way that for e = 1, we solve the 
nonlinear Schrodinger equation, and in the semiclassical limit e — > 0, we retrieve 
the limit in terms of compressible Euler equation, as recalled below. This equation 
appears in several contexts in Physics. For instance, in the case e = 1, (II. lj) 
corresponds to an envelope equation in the propagation of lasers, a case where t 
does not correspond to time, but to the direction of propagation; see e.g. [33] and 
references therein. The semiclassical regime is present in the modeling of Bose- 
Einstein condensation, where e corresponds to the (rescaled) Planck constant; see 
e.g. |40j and references therein. A remarkable property in the semiclassical regime 
is that the limit is expressed in terms of a compressible, isentropic Euler equation. 

A popular way to relate the semiclassical limit to fluid dynamics is the use of 
the Madelung transform [32 , which is essentially the polar decomposition: seek the 
solution to (jl.lj) of the form 

u e (t, x) = ^p(t,x)e iS ^' x ^ e , p^O, S e R. 
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Plugging this expression into (jl.lj) . and separating real and imaginary parts yields 
(1.2) 



y 7,[d t S+^S\ 2 + p\ =yA(Vp), 



Two comments are in order at this stage: the first equation shows that S depends on 
e and the second equation shows that so does p in general. We shall underscore this 
fact by using the notation (S £ ,p £ ). Second, the equation for S £ can be simplified, 
provided that p e has no zero. Introducing the velocity v £ — VS* 6 , (|1.2j) yields the 
system of quantum hydrodynamics (QHD), see also [21] : 



(1.3) 



d t v £ + v £ ■ Vv £ + Vp £ 
d t p £ + div (p £ v £ ) = 0. 




The term on the right hand side of the equation for v £ is classically referred to 
as quantum pressure. In the limit e — > 0, this term disappears, and we find the 
compressible Euler equation: 



(1.4) 



d t v + v ■ Vv + Vp = 0, 
d t p + div (pv) = 0. 



This approach was used recently to develop an asymptotic preserving scheme for 
the linear Schrodinger equation (lu^pu 6 is replaced with V{x)u £ ), see [17] . The goal 
of an asymptotic preserving scheme is to have a unified way to compute the solution 
as e = 1, and to retrieve the limit as e — > 0, in such a way that the discretization 
does not depend on e; see e.g. (28l [18]. As pointed out in [17], the drawback of 
Madelung transform is that it does not support the presence of vacuum (p = 0). 
The point of view that we shall study numerically is due to E. Grenier [27], and 
consists in seeking u £ as 

(1.5) u £ (t,x) = a £ (t,x)e lc/> ^ t ' x)/£ , a £ 6 C, (f> £ 6 R. 

Allowing the amplitude a £ to be complex-valued introduces an extra degree of 
freedom, compared to the Madelung transform. The choice of Grenier consists in 
imposing 

d 4 «f + i|v,f| 2 + K| 2 = o, 



(1.6) 

In terms of v £ = V0 e , this becomes 
(1.7) 



d t a £ + Vcff ■ Va £ + ^a £ A0 e = i^Aa £ 




In this model, the presence of vacuum (a £ — 0) is not a problem. We will see that 
this is so both on a theoretical level and in computational tests. In the limit e — > 0, 
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we find formally 

Vial 2 = 0, 




(L8) < " ■ " 1 rV 

-adivt) = (J. 
2 

We check that (p, v) = (|a| 2 , v) then solves (|1.4I) : (|1.8I) corresponds to the nonlinear 
symmetrization of ((0)1 f pfl[T4] ). 

In this paper, we have chosen to focus on the defocusing cubic nonlinearity, for 
which the relevance of (|1.6|) to study the semiclassical limit is proved (see H2.ip . 
It seems very likely that equivalent numerical results should be available for other 
nonlinearities, as discussed in $'2Al even though in several cases, no theoretical 
result is available concerning the natural generalization (|2.6I) of (|1.7I) . Similarly, in 
the linear setting considered in [17] , this modified Madelung transformation should 
overcome the problem of vacuum pointed out in |17) . 

We also stress the fact that the convergence of (|1.7|) towards (|1.8|) holds so long as 
no singularity has appeared in the solution of (|1.8|) (or, equivalently, in (11.41) ), Note 
that except in the very specific case d — 1 (where the cubic Schrodinger equation 
is completely integrable), no analytical result seems to be available concerning the 
asymptotic behavior of u £ as e — > for large time (that is, after a singularity has 
formed in the solution to the Euler equation) . As pointed out in [TT] , the notion of 
caustic seems to be different in the case of (|1.1[) . compared to the linear case 

e 2 

where several computational results are available past caustics (see e.g. [3SJ[25] and 
references therein, and ^2.2[) . 

1.1. Conserved quantities. Equation (|1.1|) enjoys a bi-Hamiltonian structure, 
and therefore has two quantities which are independent of time: 

d 
d 

— y\\tvu W||£2(Rd) -r ||'u ^M\L i (R^) j 

A third important quantity is conserved, which plays a crucial role, e.g. in the 
study of finite time blow-up in the case of focusing nonlinearities: 

d f 

(1.11) Momentum: — Im / u £ (t, x)eVu £ (t, x)dx = 0. 

dt J R d 

Plugging the phase/ amplitude representation (|1.5|) into these conservation laws, 
and passing formally to the limit e — > 0, we recover conservation laws associated to 
the Euler equation {HI ([TT)): 

4/ p(t,x)dx = ^- [ (p\v\ 2 + p 2 ) {t,x)dx = / (pv) (t, x)dx = 0. 
at J j^d at J p^d at J j^d 

Setting J £ (t) = x + ietV, two other evolution laws are available: 
Pseudo-conformal: ^|| J e (t)u e \\ 2 L2{Tld) + t 2 \\u £ \\ 4 Li{Ild ^ = t(2 - d^u 6 ^^. 

^-Re / u £ {t,x)J e (t)u £ {t,x)dx = 0. 
at j R d 



(1.9) Mass: -|| U £ (i)||| 2(Rti) = 0. 

(1.10) Energy: ^ f |kV U £ (t)|| 2 2(R<J) + ll^Wlli*^^ = 0. 
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Passing formally to the limit e — > 0, we infer: 

^- [ (\x-tv{t,x)\ 2 p(t,x) +t 2 p 2 (t,x)) dx = (2-d)t [ p 2 (t,x)dx. 

at J R d v / j R d 




(x — tv(t, x)) p(t, x)dx = 0. 



We discuss this aspect further into details in §2.31 

1.2. Semiclassical limit for NLS: numerical approach. The most reliable ap- 
proach so far to study numerically the semiclassical limit for Schrodinger equations 
seems to be the time-splitting spectral discretization (Lie or Strang splitting, see 
[8]): one solves alternatively two linear equations, 

e 2 

ied t v £ + — Av £ = 0, and ied t v £ = \v E \ 2 v e . 

Despite the appearance, the second equation is linear, since in view of the gauge 
invariance, d t (\v e \ 2 ) — 0, so the second equation boils down to id t v £ = \vf\ 2 v £ , 
where v\ denotes the initial value for v £ . 

Note that from [34], usual finite-difference schemes for the linear Schrodinger 
equation may lead to very wrong approximations. Instead, schemes based on the 
fast Fourier transform (FFT) have been preferred. In [5], it was shown that the 
time-splitting method, coupled with a trigonometric spectral approximation of the 
spatial derivative, conserves the total mass, and is gauge- invariant, time-reversible. 
Moreover, with this approach, the convergence of the scheme in L 2 is proved, when 
the nonlinearity in Ijl.ljl is replaced by an external potential. This regime turns 
out to be far less singular in the limit e — > than the nonlinear case of (II. ip . as 
discussed below. 

We briefly point out that the numerical study in [5J [5J shows that, contrary 
to the case of the linear Schrodinger equation, to study the semiclassical limit 
for (jl.ip with time-splitting, it is necessary to consider mesh sizes and time steps 
which are 0(e). This is due to the fact that the semiclassical regime is strongly 
nonlinear (supercritical, in the terminology of [11]): we consider initial data which 
are 0(1) in L 2 n L°°, and there is no power of e in front of the nonlinearity. As a 
consequence, the semiclassical limit is a "strongly nonlinear" process, since starting 
with a semilinear Schrodinger equation (for fixed s > 0), we come up in the limit 
e — > with a quasilinear equation (the compressible Euler equation). 

In [B], it is shown that mesh sizes and time steps must be taken of order 0(e), 
even to recover the behavior of two physically important quantities: 

Position density: p E {t,x) = |M e (i,x)| 2 = \a e {t,x)\ 2 . 
Current density: J £ (t,x) — elm(u £ (t,x)Vu £ (t,x)) . 

We refer to the numerical results in |6j Example 4.3], which show some important 
instability in the numerical approximation for (jl.ip . at least if the time step is large 
compared to e: evidently, the position and current densities cannot be computed 
correctly if mesh size and time step are independent of e. 

On the contrary, we obtain a good description of p £ and J £ as e — > when 
studying numerically the system (|1.7|) . even if the time step is independent of e. 
Things would probably be similar in the case of the QHD system (|1.3j) . up to the 
important aspect that the presence of vacuum (p £ — 0) is not allowed in (|1.3|) . 
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The idea to explain this difference is the following. To construct directly the wave 
function u £ solving errors which are large compared to e (say of order e a , 

< a < 1) lead to instability of order 0(1) on u £ after a short time (of order 
e l ~ a ). Among possible sources of errors, we can mention a simple space shift, 
which is rather likely to occur in numerical studies. This can actually be proved 
thanks to the approach of Grenier, see .10 . This is due to the strong coupling 
phase/amplitude in (|1.6j) : a small modification of the amplitude a 6 leads to a 
modification of the same order for cjf . To recover u £ , one has to divide (jf by e, 
which is small, so the actual error for <j) e may be dramatically increased. 

One can rephrase the above analysis as follows. The semiclassical limit for (11.11) is 
"strongly nonlinear" : as e — > 0, we pass from a semilinear equation (for fixed e, the 
Cauchy problem for (ll.l[) is handled by perturbative methods relying on properties 
of the linear equation, see e.g. [13]), to a quasilinear one, the Euler equation 
(jl.4[) (in which the nonlinear terms cannot be treated by perturbative methods, 
see e.g. [44] ) . As a consequence, the asymptotic behavior of u e is very sensitive 
to small errors [10] . In time splitting methods, one considers the nonlinearity as 
a perturbation, while this is not sensible in the framework of (ll.ip . unless a high 
precision in the space and time steps is demanded. It would be quite different with 
some positive power — at least 1 — of e in front of the nonlinearity; see |11) for 
theoretical explanations, and [6] for numerical illustrations. 

If one is interested only in the position and current densities, small errors in (jl.7[) 
are not so important, since one never has to divide the phase by e (see Section [2] for 
more details). This explains why we can obtain satisfactory results by considering a 
mesh size h = Ax independent of e, and a time step given by the parabolic scaling, 
that is, proportional to h 2 . 

An extra step in the numerical analysis of nonlinear Schrodinger equations was 
achieved in [7], where a semi-discrete scheme was introduced, which turns NLS 
into an almost linear system, in the case e = 1. It is based on a central-difference 
approximation shifted by a half time-step. For t n = nSt and i„+i/2 = (n + ^)5t, 
let m" be the approximation at t — t n . The scheme is given by 

/2 ( u^+u" 



n A 7T = V> 



5t 2 V 2 J r V 2 

lu"l 2 . 



^n+1/2 _j_ _0n-l/2 



This approach has the advantage of preserving the mass (|1.9j) and an analogue of 
the energy (|1.10|) of the solution [7] : 



\u n \ 2 = \u°\ 2 , and E n = E°, where 
i: •/ n. ' 

E n = J (jV-a n | 2 + 2\u n \ 2 ^ n ~ 1 ' 2 ~ (V"- 1/2 ) 



It does not seem that there is also an analogue of the momentum which is conserved, 
in the same fashion as (|1.11[) . Note that to adapt this approach numerically in the 
semiclassical regime, one would also have to consider mesh sizes and time steps 
which are O(e). Therefore, the approaches in El [7] do not seem well suited for 
asymptotic preserving schemes. 
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In this paper, we present numerical experiments only, and do not claim to justify 
the approach by numerical analysis arguments. In view of the little knowledge that 
we have on the behavior of the solution to (|1.7[) past the critical time for the Eulcr 
equation, such a study could reasonably be expected only so long as the solution 
of the Euler equation remains smooth. Yet, such a study would be an interesting 
challenge, which we do not address here. 

1.3. Outline of the paper. In Section [21 we recall the main theoretical results 
established for the semiclassical analysis of (|1.1|) . The main goal is to state some 
results which can thereafter be tested numerically to validate the scheme. The 
numerical implementation is presented in Section [3] Numerical experiments (based 
on three examples) are discussed in Sectional We conclude the paper in Section [5] 

2. The theoretical point of view 
We will always consider initial data of the form 

(2.1) u £ (0, x) = a Q {x)e l ' l ' o{x)/e , a Q e C, <£ € R, 

where ao and cf>o are smooth, say in H s (Il d ) for all s. In that case, (|1.6p is supple- 
mented with the Cauchy data 

a e (Q,x) = a (x) ; (j> £ (0, x) = <j>o(x). 

This implies that the Cauchy data for (II . T[) are 

(2.2) a e (0,x) = a (x) ; v £ (0, x) = V^,(a:). 

2.1. Known results. A second advantage of the system (|1.7|) over (|1.3|) . besides 
the role of vacuum, is that it already has the form of an hyperbolic symmetric 



system. Separate real and imaginary parts of a £ , a 6 = a\- 



(|1.7[) takes the form 



with u e 



/ of \ 



L = 






-A 





A 

















Jdxd 



V «2 / 

d 

and A(u,0 = JlA j (u)£ j = 



3=1 











«l it 
2 S 



2ai ^ 2a 2 £, v ■ £I d 

The matrix A is symmetrized by a constant diagonal matrix S such that SL = L. 
We note that L is skew-symmetric, so it plays no role in energy estimates in Sobolev 
spaces H s (R d ). The main results in [27] can be summarized as follows: 

Theorem 2.1 (From [27]). Let d ^ 1, s > 4 + d/2, and a , V0 O € H s (R d ). 

1. There exist T > tmd a unique solution (p, v) G C([0, T]; i? s (R d )) 2 to (II -4[) suc/i 
i/iai p(0,x) = |ao(cc)| 2 and v(0,x) = V<f)Q(x). 

2. For the same T, |T7]) has a unique solution (a E ,v e ) £ C([0, T]; H s ~ 2 (R d )) 2 such 
that a e (0, x) — ao(x) and v e (0, x) — \?(/)o(x). 

2'. For the same T, (| 1 . 8[) has a unique solution (a,v) G C([0,T}; H s (R, d )) 2 such 
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that a(0, x) = ao(x) and v(0, x) = V</>o(x). 
3. As e — > 0, we have: 

\\a £ - a|U°o( [0 ,T];Hs- = ) + \\V £ - u|Uoo([0,T];ff— = ^ 0?) • 

Remark 2.2 (Periodic case). The same result holds in the periodic setting (x € T d 
instead of x G R rf ), with exactly the same proof. 

Once v £ is constructed, there are at least two ways to get back to (jf . Either 
argue that v £ remains irrotational, or simply define (jf as 

(2.3) (j> £ (t) = O - j Qk e ( S )| 2 + V|a £ ( S )| 2 ) ds, 

and check that d t {v £ - V0 £ ) = d t v £ - Vd t (j> £ = 0. So for t € [0, T], that is so long 
as the solution to the Euler equation (|1.4p remains smooth, the solution to 
with initial data ttu_ = aoe^°/ £ is given by u e = a £ e l ^^ £ . 

Note that even if 0o = 0, (jf (as well as v £ and u) must not be expected to be 
zero (nor even small), because of the strong coupling in (|1.6p . Typically, if 4>o = 0, 
(H^D yields d t tf t=0 = -KIV 0. 

In addition, in the limit e — > 0, we recover the main two quadratic observables: 

p £ = \a £ \ 2 — > |a| 2 = p in L°°([0, T]; L 1 (R rf )), 

J £ = Im {eu £ Vu £ ) = \a £ \ 2 v £ + elm(a £ Wa £ ) ~~>\a\ 2 v = Jin £°°([0, T]; L x (R d )). 
We have more precisely: 

( 2 -4) \\P 6 - P||z,=°([0,T];L I n£ o °) + — ^||i oo ([0,T];Z, 1 ni«') = C (s) • 

We can also prove the convergence of the wave function (|27|). In the particular 
case which we consider where the initial amplitude a e (0,x) does not depend on e, 
we have (with an obvious definition for </>): 

\\u £ -ae^ £ \\ L ^ mnL 2 nLao) =<D(e). 

In general, a modulation of a must be taken into account to have such an approx- 
imation of the wave function ([11]): u £ ~ ae 1 ^ ] e l ^l £ . In the framework of this 
paper, we have (jy- 1 ' = (see |11[ Section 4.2]). 

2.2. An open question. As pointed out in the introduction, no analytical result 
seems to be available concerning the semiclassical limit of (|1.1[) when the solution of 
the Euler equation (|1.4p has become singular. Theorem l2.ll gives a rather complete 
picture for the asymptotic behavior of u £ for t £ [0,T], that it before the solution to 
(|1.4|) becomes singular. Note that if for instance ao and 4>q are compactly supported, 
then no matter how small they are, (a, v) develops a singularity in finite time 
( |331 IT41 141)]). On the other hand, for fixed e > 0, we know that the solution to 
(jl.ljl with initial data uf t=0 = aoe^ ^ <E H s (R d ), s ^ 1, is global in time with the 
same regularity, at least if d 4: u £ € C([0, oo[; /F(R d )). See [24] (or [13]) for the 
case d < 3, and [41] for the case d = 4 (which is energy-critical). 
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A natural question is then: what happens to u s as the solution to the Euler 
equation (II. 4|) becomes singular? In the linear setting, 



2 



(2.5) ied t uf in + -At&, = ; < in|t=0 = a e^, 

the question is rather well understood: when the solution to the corresponding 
Burger's equation (for the phase) becomes singular, a caustic is formed, which is 
a set in (i, x)-space (see e.g. [TH|3S]). Near the caustic, the amplitude of uf in is 
amplified, like a negative power of e. For instance, if 4>o(x) — —\x\ 2 /2, then 



'lin 



{t,x) 



I an ( JOJ) A*\ 2 /{Mt-i)) if + < i 




where ao denotes the Fourier transform of ao; see [TT] for several developments 
around this example, and [12] for corresponding numerical experiments. Such a 
concentration is ruled out in the case of , since the conservation of the energy 
(jl.lOp yields the uniform bound 

||it e (t)||i4( R d) ^ C independent of e (=]0, 1] and (eR. 

We remark that multiplying each equation in (jl.7[) . derivatives become exactly 
e-derivatives: every time a term is differentiated, it is multiplied by e. This is 
consistent with the possibility that a E and v E become oscillatory past the critical 
time for the Euler equation (with wavelength of order e or more). The numerical 
experiments we present below suggest that this is indeed the case. We insist on the 
fact that no result is available, though, on global existence aspects for (|1.7|) : the 
solution may be globally smooth (and e-oscillatory, in the sense of [23]), but it may 
blow up in finite time. 

Note however that the approach we present here is no longer expected to be 
asymptotic preserving beyond the breakup time for the Euler equation. The pres- 
ence of rapid oscillations is a possible explanation, and we then recover the problem 
pointed out in [5] for pre-breakup times: rapid oscillations can be resolved only if 
time step and mesh sizes are comparable to the (small) wavelength of the wave. 
Finally, we point out that even in the linear case (see e.g. the above example), one 
cannot expect an asymptotic preserving approach to solve (|2.5p after a caustic has 
formed: near the caustic, small spatial scales must be taken into account. In the 
above example, the wave function is concentrated at scale e. A possibility to get 
an aymptotically preserving approach in the linear case would be the use of La- 
grangian integrals [19]; see [IT] for an extension in a very specific nonlinear setting. 
Note however that the definition of the Lagrangian integral depends on the initial 
phase, so this approach is more delicate to implement numerically. The if -branch 
approach would lead to similar requirements; see [91 [231 |2"B] . 

2.3. Conserved quantities. In the one-dimensional case d = 1, the cubic nonlin- 
ear Schrodinger equation (|l.ip has infinitely many conserved quantities 48 (it is 
completely integrable, see [47]). We shall not emphasize this particular case in this 
paper, and rather consider the case of a cubic nonlinearity in arbitrary dimension. 
Numerical experiments are presented in the two-dimensional case d = 2. 
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In this general case, we retain the three standard conservations: mass (|1.9[) . 
momentum (|1.11[> . and energy (jl.lOp . Writing the solution to (jl.ll) as u £ — a £ e 1 ^ l £ , 
we infer three corresponding conversation laws for the solution to (|1.7[) : 

Proposition 2.3. Let d ^ 1 and (a £ ,v £ ) € C ([0, T]; H 1 n L°°{R d )) 2 solve (ITJl) . 
TTie following three quantities do not depend on time: 

(1) TTie L?-norm of a £ : ^||a e (i)|||2 (Rti) =0. 

d f 

(2) The momentum: — / (\a £ (t, x)\ 2 v £ (t, x) + e Im (a £ (t, x)\7a £ (t, x))) dx = 0. 

(3) The energy: if v? t _ Q is irrotational, V A uf i=0 = 0, £/ien 

4 / (\eVa £ {t,x) + ia £ (t,x)v £ (t,x)\ 2 + \a e (t,x)\ 4 ) dx = 0. 

Sketch of proof. This result can be proved by using the standard regularizing pro- 
cedure and suitable multipliers. We shall just indicate the formal procedure. 

The conservation of mass is proved by multiplying the second equation in (|1.7j) 
by a £ , integrating in space, and taking the real value. 

The conservation of the momentum is obtained as follows. Multiply the equation 
for v £ by ^|a e | 2 , and integrate in space. Multiply the equation for a £ by ie\7a £ + 
a £ v £ , integrate in space and consider the real value. Summing these two relations 
yields the conservation of the momentum. 

For the energy, the procedure is similar. Note that 

d t v £ = -V 0^- + |a £ | 2 ^ , hence d t (V A v £ ) = 0. 

Therefore, if V A v? t _ = 0, then we can find <j> £ such that (<j) £ ,a £ ) solves (11.61) . 
Multiply the equation in (ff by —^dt\a £ \ 2 , the equation for a £ by iedia £ + a £ d t <fi £ . 
Sum up the two equations, integrate in space, and take the real part. □ 

2.4. About other nonlinearities. Equation (ll.lj) is the defocusing cubic nonlin- 
ear Schrodinger equation. Other nonlinearities are physically relevant too: focusing 
or defocusing nonlinearities are considered, as well as other powers, in the context 
of laser Physics (see e.g. [33]) or in the context of Bose-Einstein Condensation (see 
e.g. [THUS]), f° r instance. 

The (short time) semiclassical limit for nonlinear Schrodinger equations has been 
studied rigorously for other nonlinearities. Typically, for defocusing nonlinearities 

ied t u £ + — Au £ = \u £ \ 2a u £ , a G N, 

a result similar to Theorem 12. II is available; see [3j [15] . However, the analysis does 
not rely on an extension of (jl.6|) where |a £ | 2 would be replaced with |a e | 2<T : for 
e = (corresponding to the limiting Euler equation in the case a € N too), one 
uses a nonlinear symmetrizer (the "good" unknown is (V0, a CT )), and for e > 0, 
this change of variable affects the skew-symmetric term ieAa £ in such a way that 
apparently the analysis of [37] cannot be directly adapted. 

For focusing nonlinearities, typically 

F 2 

ied t u £ + —Au £ = -\u £ \ 2a u £ , a € N, 
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the limiting equation in the system analogous to (|1.4|) is elliptic (as opposed the 
hyperbolic system (|1.4j) ). It turns out that in this case, the "elliptic Euler sys- 
tem" is ill-posed in Sobolev spaces f|37j): working with analytic regularity becomes 
necessary [37] , and sufficient [22] [45] in order to justify the semiclassical analysis. 

An hybrid nonlinearity (neither focusing, nor defocusing) also plays a role in 
physical models: the cubic-quintic nonlinearity, 



with A £ R possibly negative. This model is mostly used as an envelope equation 
in optics, is also considered in BEC for alkalimetal gases (see e.g. [20] [H [38] ) , in 
which case A < 0. The cubic term corresponds to a negative scattering length, 
and the quintic term to a repulsive three-body elastic interaction. Justifying the 
semiclassical analysis was achieved in [3] by a slight modification of the approach 
of [27J (in a different functional framework). 

To rephrase the above discussion, the approach in to study the semiclassical 
limit for 



relies on the assumption /' > 0. However, the analysis has been carried out in 
several other situations, without considering the natural generalization of (11.71) . 



It seems reasonable to believe that even though no rigorous study for this system is 
available in general (for e > 0), this system can be used for numerical simulations. 

Finally, the approach of [27] was generalized to the case where an external po- 
tential is introduced (which may model a confining trap in the framework of Bose- 
Einstein Condensation), see [XT] , and to the case of Schrodinger-Poisson system 



One expects the oscillatory nature of the solutions to be difficult to capture 
numerically. We would like to use a stable numerical scheme with the time step 
independent of e but function of h. The scheme solves system p. 61) on coarser 
meshes than what necessary to capture all wavelengths. Therefore, the solution 
has inevitably error in it. Still we give a great deal of effort on conservation issues 
for the density, energy and momentum. The time step and mesh size being both 
independent of e, one can tackle very small e values and the scheme also works for 
e = 0(1). Obviously, if all scales are aimed at being captured, then the space grid 
size need be of order of e or less and so the time step. 

Present results show that the scheme being conservative and stable macroscopic 
quantities remain observable even when the spatial-temporal oscillations are not 
fully resolved numerically because the mesh is not enough fine to capture wave- 
lengths below h. Typically, one uses h = 0.01 for all e in a square domain of side 
one. 

In our approach, conservation is ensured by projection steps to guarantee a 
correct behavior for total position, energy and current (momentum) densities. The 





(2.6) 




m uni eh Eg - 



3. Numerical implementation 
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aim is also to show that basic numerical methods [42, 39 can be used which permits 
the adaptation of generic PDE solvers. We also point out that we have privileged 
projections which are rather cheap computationally, since they are obtained by a 
simple rescaling. 

The implementation has been done in two dimensions in space but extension to 
third dimension does not appear being a difficulty Periodic boundary conditions 
and initial data with compact support have been considered. 

Let us start with system (|1.7[) which we rewrite as: 

(3.1) d t U £ + F{U e ) = 0, 

U £ {x e fl,t = 0) = U (x), U £ (dn,t) = periodic. 

where U e = (a e ,v e ). £l/j is a discrete two dimensional square domain of side L. 
Uo(x) is a regular initial condition. For all the simulations presented in this paper 
we consider 

U a (x) = {a (x),af(x),ag(x)Y, 

with ao(x) a complex function independent of e with compact support, a is real and 
/ and g are real functions with compact support. The initial pattern is therefore 
periodic of period L in both space directions. Together with the periodicity, oscil- 
lations in space can be introduced through do, / and g. Below we show numerical 
results with two values of a. 

We consider second order finite difference discretizations of partial differential 
operators. But, the periodic boundary conditions permit the implementation of 
high order spatial discretizations as well as spectral methods. One notices that 
despite the presence of first order space derivatives, no numerical viscosity is nec- 
essary to stabilize the system both in the hydrodynamic limit and for e ^ 0. We 
therefore keep the numerical viscosity to zero for all simulations which means no 
upwinding has been used. This leads to a consistent scheme with truncation error 
in h 2 : 

F(U £ )=F h {U e )+0(h 2 ). 
We consider a simple first order explicit time integration scheme: 

(3-2) \(Uln+i/2 - UU + MUU = 0, 

Uh,o = U (x h ), U E hn+1/2 {dn h ) = periodic. 

n+ 1/2 denotes an intermediate state, before projection, where conservation is not 
guaranteed for mass and momentum. It is interesting that the approach appears 
stable even for explicit time integration. With a first order scheme in time, and 
a time step in h 2 , the time integration error will be comparable to the truncation 
error in space. 

Once n+1 / 2 is computed, one needs to project it over the admissible space to 
get n+1 based on enforcing mass, energy and momentum conservation constraints 
(see Proposition 12.31) : 

. e « _ h,2M U h,n+l/2) _ 
J l,2,3\ U h,n+l/2> U h,0) - j (rre \ ~ X ' 
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where 

-M^,n+l/2) = / \ a h,n+l/2\ 2 dx, 



h{Uf in+l / 2 )= / (l a ln+l/2| 4 + £ ^ ha £ hn+l / 2 + i a h,n+l/2 V h,n+l/2 
JQh, \ 

H U h,n+l/2) = / (K,™+l/ 2 | 2l ln+l/2 + £ Im (a^n+1/2 V <,n+l/ 2 ) ) 



dx, 
dx. 



J3 is a vector of the size d of the space dimension. This problem is overdetermined 
with essentially two variables (a e and v £ ). This overdetermination is maybe one 
reason why no numerical scheme is available for these equations verifying all con- 
servation constraints. With a e complex, there are as many variables as constraints. 
Still we did not manage to enforce at the same time the mass J\ and energy J 2 
constraints. We have chosen here to enforce J\ and J3. 
J 1 can be easily enforced in a| n+1 by simply defining: 

1/2 



u h,n+l — a h,n+l/2 I j (jje \ I 

\ 11 y U h,n+l/2> ) 

The projection aims at looking for a particular equilibrium for the constraints 
after a splitting of the variables. The above scaling suggests an a prion but natural 
splitting of the variables to be modified by each constraint. More precisely, J\ 
defines the corrections for a e h n+1 / 2 an( i the vector J3 the ones for the components 
( v h n+i/2)i °f ^ ne ve l° c ity through: 

(3.3) (I 3 y = / (|a|,„ + i| 2 «,„ + i /2 )j ■ +£lm(a|, n+ i^a^ n+1 )) dx, j = l,...,d. 
Jn h K J 

Because we are looking for a cheap projection based on scaling, we adopt the 
following corrections for each component of v e h Tl+1 y 2 : 

is i 1 e \ I ^3(^/1.0) I . 1 7 

K,n+i)j = («fc,„+i/a)j I — Jf— y 3 = l,...,d. 

Through the numerical examples below we see that these scalings are efficient in 
conserving mass and current densities. 

4. Numerical experiments 

We show the application of our projection schemes for several initial conditions. 
In the first case the current density is nearly zero and not in the second. A third case 
shows the robustness of the approach with initial vanishing a £ . We show the impact 
of the projection on the conservation of mass, energy and momentum through J\, 
J2 and J3. We will see that mass and energy cannot be both conserved at the same 
time. 

4.1. Nearly zero initial current. We consider L = 0.5, ao(x) = exp(— 80((xi — 
L/2) 2 + O2 - L/2) 2 ))(l + i) and a = 1CT 10 (hence vl Q « 0). 

Figures Q] shows the initial position and current densities. Figures [2] and H show 
the solutions at T — 0.1 sec for e = 0,0.001,0.01 and 0.1 without and with the 
projection steps. 



h(Ui Q ) 
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Figure 1. Initial position (left) and norm of the current density 
vector (right) with a = 10~ 10 . 

Figures [3] and [5] show the evolution of the constraints with time for different 
values of e without and with the projection steps. The original scheme can be 
seen being not conservative and dissipative. Of course, less dissipative numerical 
schemes could be used, but this does not remove the necessity for the projection 
step. Relative momentum constraint values appear being large, but one should 
keep in mind that these are in fact very close to zero. What is most important is 
that mass and energy constraints cannot be satisfied at the same time. This can 
also be seen in the next case with initial current density. 

An interesting indicator for the behavior of the solver is by checking if the fol- 
lowing quantity is linear in e at a given time T independent of e (see (|2.4|) ): 

(4- 1 ) \\Ph,T ~ PktWl 1 ^) + \\Jh,t - J h,T\\mn h )- 

This is shown in Figure [5] at T = 0.1 sec. The slope grows with time. 

In the same way, Figure [7] shows the dependency with respect to e for the fol- 
lowing quantity (see Theorem 12. 

(4-2) |K jT - a h , T \\mn h ) + IK,t ~ v hM\^(n h )- 

Again, the dependency is linear for small e at T = 0.1 sec. 

4.2. Non zero initial current. This is the same case as before but with a = 10~ 2 
and 

f f(x) = cxp(-80((x 1 - L/2) 2 + (x 2 - L/2f)) sin^O^), 
' \ g(x) = exp(-80((afi - L/2) 2 + (x 2 - L/2) 2 )) cos(lOxi). 

Figure [5] shows the initial position and current densities. FigureOshows the solution 
at T = 0.1 sec for e = 0,0.001,0.01 and 0.1. Figure [101 shows the evolution of the 
position density, energy and current density constraints with time for different 
values of e when only mass through I\ and the current density through vector I3 
have been maintained. 
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Figure 2. Without projection. Position (left column) and norm 
of the current density vector (right column) at T = 0.1 sec for 
(resp. from the top) e = 0, 0.001, 0.01 and 0.1 with a = lO" 10 . 
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0.2 1 ' ' 1 ' ' ' ' ' ' 1 
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0.2 - 




-1500 I 1 1 1 1 1 1 1 1 1 1 

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 

Figure 3. Without projection. Evolution in time (sec) of (resp. 
from the top) the constraints on the position density, energy and 
sum of both components of the current density for s = 0, 0.001, 0.01 
and 0.1 for an initial condition with a = 10~ 10 . 
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Figure 4. Position (left column) and norm of the current density 
vector (right column) at T = 0.1 sec for (rcsp. from the top) e = 
0,0.001,0.01 and 0.1 with a = lO" 10 . 
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400 r 



0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 



Figure 5. Without projection. Evolution in time (sec) of (resp. 
from the top) the constraints on the position density, energy and 
sum of both components of the current density for e = 0, 0.001, 0.01 
and 0.1 for an initial condition with a = 10~ 10 . For the energy, 
deviation increases with e. Larger oscillations appear with higher 
e. Large values are due to the fact that initial value of the current 
is nearly zero. 
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FIGURE 6. Linear dependency of (|4.1[) at T = 0.1 sec with respect 
to e. 

10 i 1 1 1 1 1 1 1 1 




0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 



Figure 7. Linear dependency of (|4.2[) at T — 0.1 sec with respect 
to e. 

Figures [Til and [T2l show that indicators (I4.ip and (I4.2[) are still linear with respect 
to e but on a shorter range close to zero. 

4.3. a 6 changing sign. To introduce a changing sign initial data for a e , we consider 
an initial condition given by a (x) = (exp(— 320((a;i — L/2) 2 + (x% — L/2) 2 )) — 
exp(— 320((£i — L/2) 2 + (x2 — L/2) 2 )))(l + This initial amplitude changes signs: 
the set where it is zero corresponds to the presence of vacuum in the hydrodynamical 
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Figure 8. Initial position (left) and norm of the current density 
vector (right) with a = 0.01. 

point of view. The initial current is as for the previous case with / and g given in 
(|4.3[) . Figure [TBI shows the initial position and current densities. Figure [T4l shows 
the solution at T = 0.05 sec for e = 0,0.001,0.01 and 0.1. Figure [15] shows the 
evolution of the position density, energy and current density constraints with time 
for different values of e when only mass through 1\ and the current density through 
vector I 3 have been maintained. 

Figures [TBI and [T71 show indicators (I4.1j) and (I4.2j) . 

To see the behavior of the approach after singularities have formed in the Euler 
equation (for e = 0), we show in Figure [TBI the solution at T = 0.15 sec: the solution 
for e = has become singular, while the solution for e > seems to remain smooth. 
In this case, the meaning of the figure for e = is unclear, since we know that the 
scheme has dealt with a singularity. On the other hand, rapid oscillations have 
appeared at least for e — 0.1. For e — 0.001, the map is not very smooth, as if some 
oscillations were not resolved. Recall however that the time step and the mesh size 
are independent of e: in the presence of rapid oscillations, this strategy has proven 
uncfficient in [5J , as recalled in £11.21 This may very well be the case in Figure [T51 

5. Conclusion 

We have presented a numerical implementation to compute the solution of the 
system (ll.7[) . which is a way to solve the nonlinear Schrodinger equation that is 
asymptotic preserving in the semiclassical limit. To reconstruct the wave function 
u e , the phase 4> e can be computed by a simple time integration, in view of (|2.3I) . 

The scheme used in this paper is explicit, and is therefore rather cheap on the 
computational level. It preserves the L 2 -norm of the solution to the nonlinear 
Schrodinger equation, and can be adapted in order to conserve the momentum 
as well, thanks to simple projections based on rescaling. On the other hand, the 
energy is not conserved. 

With mesh sizes and time steps which are independent of the Planck constant 
e, we retrieve moreover the main two quadratic observables (position and current 
densities) in the semiclassical limit e 0, and before singularities are formed in the 
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Figure 9. Position (left column) and current density (right col- 
umn) at T = 0.1 sec for (rcsp. from the top) e = 0, 0.001, 0.01 and 
0.1 with a = 0.01. 



limiting Euler equation, up to an error of order 0(e), as predicted by theoretical 
results. The presence of vacuum (zeroes of the position density) is not a problem 



NLS AND SEMICLASSICAL LIMIT 



21 




0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 



Figure 10. Evolution in time (sec) of the constraints on the 
position density, energy and sum of both components of the current 
density (resp. from the top) for e — 0,0.001,0.01 and 0.1 for an 
initial condition with a — 0.01. 
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FIGURE 11. Linear dependency of (|4.1|) at T = 0.1 sec with re- 
spect to e. 
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1 1 1 ' ' ' ' ' 1 

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 

FIGURE 12. Linear dependency of (|4.2[) at T = 0.1 sec with re- 
spect to e. 

in this approach; the case treated in Section 14.31 is in perfect agreement with this 
theoretical result. 

Finally, these experiments suggest that once the solution to the Euler equation 
has developped singularities, the solution to (|1.7|) may remain smooth, while it 
becomes rapidly oscillatory. It is possibly e-oscillatory in the sense of [23] , but the 
existence of intermediary scales of oscillation cannot be a priori ruled out. We do 
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Figure 13. Initial position (left) and norm of the current density 
vector (right) with varying sign initial a £ . 

not claim to observe any quantitative result for post-breakup time, but rather a 
qualitative phenomenon: a refinement of time step and mesh size would be needed 
in view of a more reliable result after the breakup time. This aspect goes beyond 
the scope of the present paper. 
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